function auto_bridge();

% auto_bridge.m : take data traces in current clamp with a pulse, and automatically
% attempt to compute the optimal bridge balance for each trace.
% Once computed, the bridge values are stored in CONTROL(sel).Bridge, and can be applied
% when the data is read back into the program

% if(~isempty(tfit) & any(strcmp(protocol, {'iv', 'hyp2', 'ivf', 'cciv', 'ap-iv', 'ap_iv', 'ap-iv2'}))) % have data to use


%try
sf = getmainselection;
spike_flag = 0;
if(sf > 0) 
   pflag = getplotflag;
   QueMessage('Auto_Bridge', 1); % clear the que
   for i = 1:length(sf)
      [dr{i} tr{i}] = autobr(sf(i), pflag, spike_flag);
   end;
end;
figure;
ts = 0;
for i = 1:length(sf)
   plot(tr{i} + ts, dr{i}, 'ko');
   ts = max(tr{i});
   hold on
end;
%catch
%  watchoff;
%   QueMessage('Error in Current-clamp Analysis routine', 1);
%end;

function [dr, tr] = autobr(sf, plot_flag, n)

global VOLTAGE CURRENT DFILE
global CONTROL

err = 0;

QueMessage('Auto_Bridge - Starting');
dat = [];
time = [];
ivresult = []; % initialize it
do_spike = 1;
if(nargin == 3)
   do_spike = number_arg(n)-1;	% 1 turns off detailed spike analysis; 2 turns it on
else
   do_spike = 0; % assume we do not want it
end;

[DFILE, err] = analysis_setup(DFILE, sf);

if(err ~= 0)
   return;
end;



% abstract general information
deadwin = CONTROL(sf).deadwin; % dead window for most analyses EXCEPT spikes, which have their own values.

spike_thresh=number_arg(CONTROL(sf).thresh);
QueMessage(sprintf('Auto_br: Spike Threshold: %7.2f mV', spike_thresh));
protocol=deblank(lower(CONTROL(sf).protocol));

[records,pts]=size(CURRENT);

% get windows 
ts1=number_arg(CONTROL(sf).durho);
ts2=number_arg(CONTROL(sf).durs1);
ts3=number_arg(CONTROL(sf).durs2);
ts4=number_arg(CONTROL(sf).durs3);

% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
RATES = (DFILE.rate .* DFILE.nr_channel) / 1000; % array of sampling rates, convert to msec

tdel = ts1; % for auto_bridge, we use the first step
tstepdur = ts2; % and duration.
% mean times:
m1=tdel+tstepdur*0.5;
m2=tdel+tstepdur*0.99;
m3=tdel+tstepdur*0.97; % use for h current assessment

%peak times:
p1=tdel+tstepdur*0.01;
p2=tdel+tstepdur*0.50;
p3=tdel+tstepdur*0.25;


% calculate measurement times
trmp = floor(ts1./RATES); % for RMP/ihold determination

mean_t1=floor(m1./RATES); % for ss iv
mean_t2=floor(m2./RATES);
mean_t3=floor(m3./RATES);

peak_t1=floor(p1./RATES); % for peak iv
peak_t2=floor(p2./RATES);
peak_tfit=floor(p3./RATES);
end_tfit = floor(m2./RATES);
pre_g2 = floor((tdel-2)./RATES); % pre membrane potential measurement...

% but first smooth the voltage out a bit
for i = 1:records
   fsamp = 1000/RATES(i); % get sampling frequency
   fco = 3000;		% cutoff frequency in Hz
   wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
   if(wco < 1) % if wco is > 1 then this is not a filter!
      [b, a] = Butter(8, wco); % butterworth
      vsmo(i,:) = filter(b, a, VOLTAGE(i,:)); % filter all the traces...
   else
      vsmo(i,:) = VOLTAGE(i,:);
   end
end

QueMessage('AutoBridge - Finding mean and peak voltage');
% measure holding current and "rmp"
for i=1:records
   hold_cur(i) = mean(CURRENT(i,1:trmp(i))); % save in array for later usage
   rmp_volt(i) = mean(VOLTAGE(i,1:trmp(i)));
   
   % measure steady-state IV and voltage
   mean_volt(i) = mean(VOLTAGE(i,mean_t1(i):mean_t2(i))');
   final_volt(i) = mean(VOLTAGE(i,mean_t3(i):mean_t2(i))');
   mean_cur(i) = mean(CURRENT(i,mean_t1(i):mean_t2(i))');
end

QueMessage('autobr - time base 2');
if(DFILE.mode >= 5)
   wz=DFILE.ztime;
   w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
   if(~isempty(w))
      wz(w+1:end)=wz(w+1:end)+24*60*60;
   end;
   zt = (wz-wz(1))/(60);
   cond_baseline = 5; % 5 min baseline
else
   zt = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
   cond_baseline = 3;
end;
tr=zt;


QueMessage('autobr: Computing taus');

for rec=1:records
   g1(rec) = floor(tdel/RATES(rec))-1;
   if(~isempty(deadwin))
      g2(rec) = floor((tdel+deadwin)/RATES(rec));
   else
      g2(rec) = g1(rec);
      deadwin = 0;
   end
end
%figure;
for rec=1:records
   this_rmp(rec) = mean(VOLTAGE(rec,pre_g2(rec):g2(rec)-1));
   i0=mean(CURRENT(rec,pre_g2(rec):g2(rec)-1));
   i1=mean(CURRENT(rec, g2(rec):end_tfit(rec)));
   t_list(rec,:) = [g2(rec):end_tfit(rec)]; 
   time_fit = time(rec, t_list(rec,:))-tdel-deadwin;
   % do a single exp fit first
   [a0 a1 tau] = expfit(time_fit, vsmo(rec, t_list(rec,:)));
   if(a1 < 0 | tau < 0)
      a1 = 5; % fake it...
      tau = 5;
   end      
  % hold off
   trace_index = [pre_g2(rec):end_tfit(rec)];
 %  plot(time(rec, trace_index)-tdel-deadwin, vsmo(rec, trace_index));
  % hold on
   %   [a0 a1 tau]
   [fpar, chisq, niter, volt_fit1(i,:)] = mrqfit('exponential', [a0 a1 tau], time_fit, vsmo(rec, t_list(rec,:)), [], ...
      [], [], [], 50, []);
   if(isempty(volt_fit1))
      volt_fit1(i,:) = a0 + a1*exp(-time_fit/tau);
   end
   fp1(i,:)=fpar;
   cerr1(i)=chisq;
   amp0 = fpar(1); amp1=fpar(2); tau1 = fpar(3);	      
   [fpar2, chisq2, niter2, volt_fit2(i,:)] = mrqfit('exponential', [amp0 amp1/2 tau1 amp1/2 tau1/8], time_fit, vsmo(rec, t_list(rec,:)), [], ...
      [], [], [], 50, []);
   amp2 = fpar2(1); amp21=fpar2(2); tau21 = fpar2(3); amp22 = fpar2(4); tau22 = fpar2(5);
   %time_fit = time_fit+tdel+deadwin;
   if(tau21 > 50 | tau22 > 50)
      chisq2=5*chisq; % make it a worse fit.
   end
   fp2(i,:)=fpar2;
   cerr2(i)=chisq2;
	estrmp = amp2+amp21+amp22;
	dr(rec) = (this_rmp(rec) - estrmp)/(i0-i1);
%	fprintf(1, 'i0: %7.3f i1: %7.3f Bridge Error: %7.3f Meghoms\n', i0, i1, -1000*dr(rec));
%   vcorr = [(VOLTAGE(rec, pre_g2(rec):g2(rec)-1) - i0*dr(rec)) (VOLTAGE(rec, g2(rec):end_tfit(rec))-i1*dr(rec))];
%      plot(time_fit, volt_fit2(i,:), 'r');
%   plot([min(time(rec, trace_index)-tdel-deadwin) max(time(rec, trace_index)-tdel-deadwin)], [estrmp, estrmp], 'k');
%   plot(time(rec, trace_index)-tdel-deadwin, vcorr, 'g');
%   drawnow
%   pause(0.1);
end % of for loop over records to be fit


return;

